HFMD 2023 outbreak in Ho Chi Minh City

Author

Nguyen Pham Nguyen The

Published

07-16-2025

1 Background

A hand-foot-and-mouth disease (HFMD) severe outbreak stroke Ho Chi Minh in 2023, with two successive peaks in incidence in July and September. Laboratory analyses showed that both peaks were caused by enterovirus A71 (EV-A71) subgenogroup B5 lineage.

We consider 3 hypotheses that could explain this pattern:

  • The peaks were caused by cases from locations

  • The 2 successive peaks reflected the spread of the virus across age classes

  • The period from July to August coincides with summer break characterized by a drop in the contact rate.

2 Data analysis

2.1 Data structure

From January to December 2023, a total of 43,372 HFMD cases were reported in Ho Chi Minh City.

Characteristic 2013
N = 8,0781		 2014
N = 10,0431		 2015
N = 8,7291		 2016
N = 5,7401		 2017
N = 30,8251		 2018
N = 39,3521		 2019
N = 28,4861		 2020
N = 16,3981		 2021
N = 9,7721		 2022
N = 19,1651		 2023
N = 43,3801		 2024
N = 7,3971
Age NA (NA, NA) NA (NA, NA) NA (NA, NA) NA (NA, NA) 1.84 (0.69, 5.06) 1.95 (0.69, 5.68) 1.85 (0.66, 5.44) 1.83 (0.71, 5.02) 2.21 (0.74, 5.45) 2.25 (0.78, 5.50) 2.54 (0.79, 6.04) 2.35 (0.69, 5.95)
Gender











    Female 3,159 (39%) 4,070 (41%) 3,617 (41%) 2,391 (42%) 12,982 (42%) 16,815 (43%) 12,196 (43%) 6,944 (42%) 4,123 (42%) 7,950 (41%) 18,461 (43%) 3,121 (42%)
    Male 4,919 (61%) 5,973 (59%) 5,112 (59%) 3,349 (58%) 17,843 (58%) 22,537 (57%) 16,290 (57%) 9,454 (58%) 5,649 (58%) 11,215 (59%) 24,919 (57%) 4,276 (58%)
1 Median (5% Centile, 95% Centile); n (%)

2.2 Hypothesis 1: The peaks were caused by cases from locations

We looked at epicurves per commune and/or district (i.e. whether some locations have the first peak only and some other locations have the second peak only)

Epidemic curve by admission date for 22 districts in Ho Chi Minh City.
Conclusion

All districts exhibited a two-peak epidemic curve, so we rejected this hypothesis.

2.3 Hypothesis 2: The 2 successive peaks reflected the spread of the virus across age classes

We looked at the age distribution of cases as a function of time and space (commune or district). In this analysis we considered 2 options to deal with time: either as a continuous variable, or as a binary variable (first peak vs second peak).

Age density was estimated using kernel density estimation with Gaussian kernel density estimator, as follows:

\[ \hat{f}_x(h) = \frac{1}{n h \sigma \sqrt{2\pi}} \sum_{i=1}^{n} \exp\left( -\frac{(x - x_i)^2}{2 h^2 \sigma^2} \right) \]

Where n is the sample size, h is the bandwidth, and σ is the standard deviation of the sample. We follow Silverman’s “rule of thumb” for choosing the bandwidth of a Gaussian kernel density estimator (Silverman 2018).

\[ h = 0.9An^{\frac{-1}{5}} \] The Wilcoxon test was used to compare age distribution between the two peaks, p-value < 0.05 was considered statistically significant.

Based on the trough between the two peaks in week 36 (September 6, 2023), we divided the outbreak into two waves: the first from January to September 6, 2023, and the second from September 13 to December 2023, with 20,358 and 23,014 HFMD reported cases, respectively.

(A) Weekly case count of HFMD outbreak in Ho Chi Minh City in 2023, the vertical line separates the cases of the first and the second wave of the outbreak for analysis in Figure 1C. (B) Weekly distribution of age over time, the lines represent the baby cohorts through a year. (C) Age distribution between the first and second wave of the outbreak.
Characteristic 1st peak
N = 20,3581		 2nd peak
N = 23,0141		 p-value2
age 2.649 (1.784, 3.726) 2.393 (1.548, 3.721) <0.001
1 Median (Q1, Q3)
2 Wilcoxon rank sum test
Results

  • The first wave was caused mainly by two age groups: 1–2 years old and 2–3 years old, while the 0–1 year old and 1–2 years old groups contributed to the second half

  • The newborn baby cohort showed the first signal of infection in September 2023 and mainly contributed to the second wave from October to December 2023

Conclusion

  • The outbreak was caused by two different age groups: children aged 1–3 years mainly contributed to the first wave, while those aged 0–2 years mainly contributed to the second wave

  • The low density of the 0–1 age group during the first peak, when these children were around six months old, indicates that maternal immunity offered protection during the first six months of life

2.4 Hypothesis 3: Variation in contact rate between summer break and school term

We used the time-series SIR (tSIR) model (Becker and Grenfell 2017) to evaluate the interaction between the level of population susceptible and contact rate.

\[ \begin{align*} S_{t+1} &= B_{t+1} + S_t - I_{t+1} \\ \mathbb{E}(I_{t+1}) &= \beta_{t+1} S_t I_t^{\alpha} \end{align*} \]

Where:

  • S: susceptible population
  • I: Infectious population
  • B: Birth
  • \(\alpha\): epidemic saturation (correction factor for switching from continuous to discrete time)
  • \(\beta\): seasonal contact rate

Main assumption:

  • The infectious period is fixed at the sampling interval of the data (weekly)

  • Over a long enough time, the sum of births and cases should be approximately equal

2.5 Birth estimation

The Hepatitis B vaccine registry of HCMC was used for birth data of the tSIR model

General Additive Model (GAM) with REML method for births estimation in 2023:

\[ Births \sim s(week_{1-52}) + s(week_{1-270}) \]

2.6 tSIR model

3 Serological data analysis

3.1 Seroepidemiology

The generalized additive models (GAM) with a binomial error distribution, logit link, and B-splines smooth terms with the restricted maximum likelihood method were used to model seroprevalence of each age group as a function of time, and seroprevalence as a function of age.

3.2 Compare attack rate and seroprevalence per age group over time

Attack rates by age group were calculated using the population size for each age group from the latest Ho Chi Minh city census in 2019. Attack rate (AR) on day t was followed as:

\[ AR_{t} = \frac{Cases_{t}}{Population \ size - \sum_{i=1}^{t-1} Cases_{i}} \]

4 Temporally Forced Model

4.1 Background

Traditional SIR model assumes that the force of infection depends on the proportion of infectives in the population, formulated as \(\beta \frac{I}{N}\) in the following system of ODEs.

\[ \begin{cases} \frac{dS}{dt} = \mu N-\beta \frac{I}{N} S\\ \frac{dI}{dt} = \beta \frac{I}{N} S - \gamma I\\ \frac{dR}{dt} = \gamma I\\ \end{cases} \]

Where:

  • \(\mu\) is the per capita birth rate

  • \(\beta\) is the transmission rate

  • \(\gamma\) is the recovery rate

Note

A range of statistical approaches have revealed that transmission of childhood infections varies seasonally, peaking at the start of the school year and declining significantly in the summer months.

4.2 Sinusoidal \(\beta(t)\)

(Bailey 1975) explored a simplified SIR model:

\[\frac{dS}{dt} = \mu N - \beta(t) \frac{I}{N} S\]

\[\frac{dI}{dt} = \beta(t) \frac{I}{N} S - \gamma I \] where:

  • \(\mu\) is the per capita birth rate

  • \(\gamma\) is the recovery rate from the infection

The transmission rate is a function of time, β(t), and was taken by Bailey to be a sinusoid:

\[\beta(t) = \beta_{0}(1+\beta_{1}cos(\omega t))\]

where:

  • \(\beta_{0}\): the baseline or average transmission rate

  • \(\omega\): the period of the forcing (\(\frac {2 \pi}{time~unit}\),1 year = \(\frac{2 \pi}{1}\))

  • \(\beta_{1}\): the amplitude of seasonality

4.2.1 SIR model with seasonality

The first systematic examination of seasonality affecting the dynamical pattern of epidemics was made, as far as we are aware, by Klaus Dietz in his seminal 1976 paper. Dietz carried out a stability analysis of the familiar SIR model:

\[ \begin{cases} \frac{dS}{dt} = bN - (\beta(t) \frac{I}{N} + b)S \\ \frac{dI}{dt} = \beta(t) \frac{I}{N} S - (b + \gamma)I \end{cases} \]

Where \(\beta(t)\) is defined as \(\beta(t) = \beta_0(1- \beta_1 \text{cos}(\omega t))\) ( he used a minus sign in his formulation in order to ensure that contact rates were at their lowest at the start of the epidemic year)

Code 
curve(1*(1 - 0.5*cos(2*pi*time)), 0, 6,
      n = 1001, xname = "time", xlab = "Time", ylab = "Transmission rate")

SIR model with seasonality 
library(odin2) 
library(dust2)

sir_seasonality <- odin2::odin({
  N <- S + I + R
  deriv(S) <- b*N - (beta_t*(I/N) + b)*S
  deriv(I) <- beta_t*(I/N)*S - (b + gamma)*I
  deriv(R) <- gamma*I - R*b

  # seasonality forcing
  beta_0 <- parameter(0.4)
  beta_1 <- parameter(0)
  omega <- parameter(2*3.14/52) # use week as time unit by default
  beta_t <- beta_0*(1 - beta_1*cos(omega*time))
  
  # initialize starting population
  init_S <- parameter(9500)
  init_I <- parameter(500)
  gamma <- parameter(0.05)
  b <- parameter(0.05)
  
  initial(S) <- init_S
  initial(I) <- init_I
  initial(R) <- 0
  
})
helper function 
# helper function to run the model 
run_mod <- function(mod, pars, duration=100, timestep=0.05){
  # --- initialize simulation time ---- 
  times <- seq(0, duration, timestep)
  
  sys <- dust_system_create(mod, pars)
  dust_system_set_state_initial(sys)
  out <- dust_system_simulate(sys, times)
  out <- dust_unpack_state(sys, out)
  
  out %>% 
    as.data.frame() %>% 
    mutate(
      t = times
    )
}

# helper function to plot the proportion of a specified compartment 
plot_comp_prop <- function(data, comp="I", t_lower = 90, use_log = TRUE, osc_freq=NULL){
  plot <- data %>% 
    mutate(
      N = S + I + R,
      prop = .data[[comp]]/N,
      prop = if(use_log) log(prop) else prop
    ) %>% 
    filter(t >= t_lower) %>%
    ggplot(aes(x = t, y = prop)) +
      geom_line(color = "cornflowerblue") 
  
  # --- annotate start and end of each cycle if oscillation frequency is provided --- 
  if(!is.null(osc_freq)){
    # compute oscillation period
    osc_period <- 2*pi/osc_freq
    
    cycle_df <- data.frame(
        cycle_t = seq(0, by = osc_period, 
                      length.out = round(max(data$t)/osc_period) + 1)
      ) %>% 
      filter(cycle_t >= t_lower)
    
    plot <- plot + geom_vline(
      aes(xintercept = cycle_t),
      data=cycle_df,
      color = "red", linetype = "dashed")
  }
  
  plot + 
    labs(
        title = paste0("Oscillations in ",comp, " compartment"),
        y = if(use_log) paste0("Log(proportion of ",comp, ")") else 
          paste0("proportion of ", comp),
        x = "Time (years)"
      )
}

# helper function to plot the SIR model
plot_sir <- function(data, ylab="Count"){
  data %>% 
    ggplot() +
      geom_line(aes(x = t, y = S, color = "S")) +
      geom_line(aes(x = t, y = I, color = "I")) +
      geom_line(aes(x = t, y = R, color = "R")) +
      scale_color_manual(
        values = c("S" = "cornflowerblue", "I" = "red", "R" = "green")
      )+
      labs(
        title = "Model output",
        x = "Time",
        y = ylab
      )
}

Model

\[ \begin{cases} \frac{dS}{dt} = bN - (\beta(t) \frac{I}{N} + b)S \\ \frac{dI}{dt} = \beta(t) \frac{I}{N} S - (b + \gamma)I \end{cases} \]

Where \(\beta(t) = \beta_0(1- \beta_1 \text{cos}(\omega t))\)

Code for SIR without seasonality forcing 
no_forcing_pars <- list(
  beta_1 = 0,
  beta_0 = 478,
  gamma = 365/13,
  b = 0.02,
  init_S = 6e-2,
  init_I = 1e-3
)

# --- run model + plot proportion of infectives ----
run_mod(sir_seasonality, no_forcing_pars, duration=150) %>% 
  plot_comp_prop(t_lower=0)+
  theme_minimal()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_text(size = 18))

Code for SIR with seasonality forcing 
forcing_pars <- no_forcing_pars
forcing_pars$beta_1 <- 0.1
forcing_pars$omega <- 3

# --- run model + plot proportion of infectives ----
run_mod(sir_seasonality, forcing_pars, duration=150) %>% 
  plot_comp_prop(t_lower=0)+
  theme_minimal()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_text(size = 18))

Conclusion

In the absence of seasonal forcing, the SIR family of models exhibit a stable equilibrium. The introduction of time-dependent transmission rates can generate a variety of dynamical patterns—depending on parameter values— ranging from simple annual epidemics to multiennial outbreaks and eventually chaos.

4.2.2 Dynamical Consequences of Seasonality

In this section, we will explore how parameter values can generate a variety of dynamical patterns

Harmonic and Sub-harmonic resonance

In the absence of seasonal forcing (i.e. \(\beta_1 = 0\)), the system fluctuates with frequency \(F\), where.

\[ F^2 = b(\gamma + b)(R_0 - 1) - (\frac{b R_0}{2})^2 \]

These oscillations are hereby referred to as the natural oscillations.

When the natural period of oscillations is approximately the same as seasonal forcing period (i.e., \(F \approx \omega\) ) we observe harmonic resonance

When natural period of oscillations \(\frac{1}{F}\) is close to an integer multiple of period of forcing \(\frac{1}{\omega}\) (i.e. \(\frac{1}{F} \approx n \frac{1}{\omega}\)), we observe a phenomenon called sub-harmonic resonance.

get_f_val function 
get_f_val <- function(params, R_0=NULL){
  f_val <- NULL
  # --- Compute F -----
  with(params, {
    R_0 <- (beta_0)/(gamma)
    f_val <<- sqrt( b*(gamma + b)*(R_0 - 1) - (b*R_0/2)**2 )
  })
  
  f_val
}
harmonic resonance 
harmonic_pars <- no_forcing_pars
harmonic_pars$beta_1 <- 0.1

# compute natural oscillation frequency  
f_val <- get_f_val(harmonic_pars)
# set omega = F
harmonic_pars$omega <- f_val

# --- run model + plot proportion of infectives ----
run_mod(sir_seasonality, harmonic_pars, duration=150) %>%
  plot_comp_prop(t_lower=140)+
  theme_minimal()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_text(size = 18))

sub-harmonic resonance 
sub_harmonic_pars <- no_forcing_pars
sub_harmonic_pars$beta_1 <- 0.1

# compute natural oscillation frequency  
f_val <- get_f_val(sub_harmonic_pars)
# set omega = F
sub_harmonic_pars$omega <- 2*f_val

# --- run model + plot proportion of infectives ----
run_mod(sir_seasonality, sub_harmonic_pars, duration=150) %>%
  plot_comp_prop(t_lower=140)+
  theme_minimal()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_text(size = 18))

Changes in \(\beta_1\) and \(R_0\) can also lead to qualitatively different epidemic patterns.

\(\beta_{1}\)
Code 
## bifurcation diagram
vlines <- c(131.5,133.2,135.65,137.40,139.85,141.60,144.05,145.80,148.25,149.95)

bete1_seq <- seq(0.1,0.5,by = 0.001)
out.df <- matrix(NA, ncol = 2, nrow = 0)
for (z in 1:length(bete1_seq)) {
  forcing_pars$beta_1 <- bete1_seq[z]
  checkdt <- run_mod(sir_seasonality, forcing_pars, duration=150) %>%
    mutate(
      N = S + I + R,
      prop = log(I/N)
    ) %>% 
    filter(t >= 70)  
  uval <- checkdt %>% filter(t %in% vlines) %>% pull(prop) %>% unique()
  out.df <- rbind(out.df, cbind(rep(bete1_seq[z], length(uval)), uval))
}
out.df <- as.data.frame(out.df)
colnames(out.df) <- c("beta_1", "N")

bifur <- ggplot(out.df, aes(x = beta_1, y = N)) + 
  geom_point(size = 1)+
  scale_x_continuous(breaks = seq(0.1,0.5,by = 0.05))+
  scale_y_continuous(limits = c(-30,0))+
  geom_vline(xintercept = c(0.1,0.32,0.37,0.45))+
  theme_classic()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_text(size = 18))

## 
forcing_pars$beta_1 <- 0.1
beta_0.15 <- run_mod(sir_seasonality, forcing_pars, duration=150) %>%
  plot_comp_prop(t_lower=130)+
  geom_vline(xintercept = vlines, color = "red",alpha = 0.2)+
  labs(subtitle = "beta1 = 0.1")+
  theme_classic()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18))

##
forcing_pars$beta_1 <- 0.32
beta_0.32 <- run_mod(sir_seasonality, forcing_pars, duration=150) %>%
  plot_comp_prop(t_lower=130)+
  geom_vline(xintercept = vlines, color = "red",alpha = 0.2)+
  labs(subtitle = "beta1 = 0.32")+
  theme_classic()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18))

##
forcing_pars$beta_1 <- 0.37
beta_0.37 <- run_mod(sir_seasonality, forcing_pars, duration=150) %>%
  plot_comp_prop(t_lower=130)+
  geom_vline(xintercept = vlines, color = "red",alpha = 0.2)+
  labs(subtitle = "beta1 = 0.37")+
  theme_classic()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18))

##
forcing_pars$beta_1 <- 0.45
beta_0.45 <- run_mod(sir_seasonality, forcing_pars, duration=150) %>%
  plot_comp_prop(t_lower=130)+
  geom_vline(xintercept = vlines, color = "red",alpha = 0.2)+
  labs(subtitle = "beta1 = 0.45")+
  theme_classic()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18))

library(patchwork)

(beta_0.15| beta_0.32 |beta_0.37|beta_0.45)/
  bifur

Constructing a bifurcation diagram. The top four panels depict time-series data for the SIR model with different levels of seasonality (beta_1 = 0.1, 0.32, 0.37 and 0.45, respectively). The vertical line at the top of the panels indicate the points when the time series are sampled in order to construct the bifurcation diagram below. The vertical line in the bottom figure indicates the value of beta_1 corresponding to the top figure.
\(R_{0}\)
Code 
## bifurcation diagram
vlines <- vlines <- c(130.95,133.70,135.15,137.90,139.35,142.05,143.5,146.25,147.7)

bete0_seq <- seq(100,1000,by = 0.5)
out.df <- matrix(NA, ncol = 2, nrow = 0)
for (z in 1:length(bete0_seq)) {
  forcing_pars$beta_0 <- bete0_seq[z]
  checkdt <- run_mod(sir_seasonality, forcing_pars, duration=150) %>%
    mutate(
      N = S + I + R,
      prop = log(I/N)
    ) %>% 
    filter(t >= 70)  
  uval <- checkdt %>% filter(t %in% vlines) %>% pull(prop) %>% unique()
  out.df <- rbind(out.df, cbind(rep(bete0_seq[z], length(uval)), uval))
}
out.df <- as.data.frame(out.df)
colnames(out.df) <- c("beta_0", "N")

bifur_r0 <- ggplot(out.df, aes(x = beta_0, y = N)) + 
  geom_point(size = 1)+
  geom_vline(xintercept = c(250,600,900))+
  theme_classic()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_text(size = 18))

##
forcing_pars$beta_0 <- 250
beta0_250 <- run_mod(sir_seasonality, forcing_pars, duration=150) %>%
  plot_comp_prop(t_lower=130)+
  geom_vline(xintercept = vlines, color = "red",alpha = 0.2)+
  labs(subtitle = "beta0 = 250")+
  theme_classic()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18))

##
forcing_pars$beta_0 <- 600
beta0_600 <- run_mod(sir_seasonality, forcing_pars, duration=150) %>%
  plot_comp_prop(t_lower=130)+
  geom_vline(xintercept = vlines, color = "red",alpha = 0.2)+
  labs(subtitle = "beta0 = 600")+
  theme_classic()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18))

##
forcing_pars$beta_0 <- 900
beta0_900 <- run_mod(sir_seasonality, forcing_pars, duration=150) %>%
  plot_comp_prop(t_lower=130)+
  geom_vline(xintercept = vlines, color = "red",alpha = 0.2)+
  labs(subtitle = "beta0 = 900")+
  theme_classic()+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18))

(beta0_250 | beta0_600 | beta0_900)/
  bifur_r0

Constructing a bifurcation diagram. The top four panels depict time-series data for the SIR model with different levels of seasonality (beta_0 = 250, 600, and 900, respectively). The vertical line at the top of the panels indicate the points when the time series are sampled in order to construct the bifurcation diagram below. The vertical line in the bottom figure indicates the value of beta_0 corresponding to the top figure.

4.3 Mechanisms of Multi-Annual Cycles

The condition for growth of disease incidence is

\[\begin{align} \frac {dI}{dt} &= \beta S \frac {I}{N} - (\gamma + \mu) I > 0 \\ &= I (\gamma + \mu) (R_{0}\frac{S}{N}-1) >0 \\ &\Rightarrow \frac{S}{N} >\frac{1}{R_{0}} \approx \frac{\gamma}{\beta} \end{align}\]

The spread only occur when a sufficient fraction of susceptible higher a critical value determined by \(R_{0}\)

Code 
harmonic_pars <- no_forcing_pars
harmonic_pars$beta_1 <- 0.1
harmonic_pars$omega <- f_val

harmonic_out <- run_mod(sir_seasonality, harmonic_pars, duration=150) %>% 
  filter(t>=140) 

harmonic_out %>% 
  mutate(
    # compute beta_t and threshold for the filtered timeframe
    N = S + I + R,
    prop_S = S/N,
    prop_I = I/N,
    beta_t = harmonic_pars$beta_0*(
      1 - harmonic_pars$beta_1*cos(harmonic_pars$omega*t)
    ),
    threshold = (harmonic_pars$gamma)/beta_t, 
    # check whether prop of S is above threshold
    above_threshold = prop_S > threshold
  ) %>% 
  filter(t >= 140) %>% 
  ggplot(aes(x = t))+
  geom_line(aes(y = prop_S),color = "cornflowerblue",lwd=1)+ 
  geom_point(aes(y = prop_S, color = above_threshold), size = 2.5) +
  geom_line(aes(y = threshold), 
            color = "black", linetype = "dashed",lwd=1)+
  geom_line(aes(y = 0.1-(log(prop_I)/-200)),lwd=1,alpha = 0.5)+
  geom_segment(aes(x = 139.5, y = 0.065, xend = 139.9, yend = 0.065),
               arrow = arrow())+
  annotate("text", x= 139.05, y=0.065, label= "Infectives",size = 5)+
  geom_segment(aes(x = 139.5, y = 0.061, xend = 139.9, yend = 0.061),
               arrow = arrow())+
  annotate("text", x= 139.05, y=0.061, label= "Threshold",size = 5)+
  geom_segment(aes(x = 139.5, y = 0.049, xend = 139.9, yend = 0.049),
               arrow = arrow())+
  annotate("text", x= 139.05, y=0.049, label= "S/N",size = 5)+
  theme_classic()+
  scale_y_continuous(name = "Susceptible")+
  scale_x_continuous(limits = c(139,150))+
  labs(color='Above threshold')+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 18))

The anatomy of seasonally forced epidemics. The level of susceptibles is color coded indicating when it is above or below the threshold.

  • The peak in disease incidence coincides with the point at which the fraction of susceptibles \(\frac{S}{N}\) falls below \(\frac{\gamma}{\beta(t)}\).

  • Once \(S > \frac{\gamma}{\beta(t)}\) disease incidence rises. This pattern is repeated once transmission has depleted susceptibles below the threshold

Code 
sub_harmonic_out <- run_mod(sir_seasonality, sub_harmonic_pars, duration=150) %>% 
  filter(t>=140) 

sub_harmonic_out %>% 
  mutate(
    # compute beta_t and threshold for the filtered timeframe
    N = S + I + R,
    prop_S = S/N,
    prop_I = I/N,
    beta_t = sub_harmonic_pars$beta_0*(
      1 - sub_harmonic_pars$beta_1*cos(sub_harmonic_pars$omega*t)
    ),
    threshold = (sub_harmonic_pars$gamma)/beta_t, 
    # check whether prop of S is above threshold
    above_threshold = prop_S > threshold
  ) %>% 
  filter(t >= 140) %>% 
  ggplot(aes(x = t))+
  geom_line(aes(y = prop_S),color = "cornflowerblue",lwd=1)+ 
  geom_point(aes(y = prop_S, color = above_threshold), size = 2.5) +
  geom_line(aes(y = threshold), 
            color = "black", linetype = "dashed",lwd=1)+
  geom_line(aes(y = 0.1-(log(prop_I)/-200)),lwd=1,alpha = 0.5)+
  geom_segment(aes(x = 139.5, y = 0.069, xend = 139.9, yend = 0.069),
               arrow = arrow())+
  annotate("text", x= 139.05, y=0.069, label= "Infectives",size = 5)+
  geom_segment(aes(x = 139.5, y = 0.055, xend = 139.9, yend = 0.055),
               arrow = arrow())+
  annotate("text", x= 139.05, y=0.055, label= "Threshold",size = 5)+
  geom_segment(aes(x = 139.5, y = 0.061, xend = 139.9, yend = 0.061),
               arrow = arrow())+
  annotate("text", x= 139.05, y=0.061, label= "S/N",size = 5)+
  theme_classic()+
  scale_y_continuous(name = "Susceptible")+
  scale_x_continuous(limits = c(139,150))+
  labs(color='Above threshold')+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 18))

The anatomy of seasonally forced epidemics. The level of susceptibles is color coded indicating when it is above or below the threshold

  • For a slightly higher amplitude of seasonality (\(\omega = 2F\)) the picture changes in important ways

  • In this instance, however, the peak in the infectives is substantially larger than in the previous graph (due to the greater transmission rate \(\omega = 2F\))

  • While at this point \(S > \frac{\gamma}{\beta(t)}\), the transmission rate is very near its annual maximum (the threshold is near its minimum), and as a result the susceptibles do not remain above the threshold long enough to produce a large epidemic. It is only when the level of susceptibles exceeds the threshold for a second time that a large epidemic begins

Conclusion

In this instance, the seasonal changes in the transmission rate are intimately associated with the dynamics of susceptibles

5 Forcing function

We’ve explored seasonality by assuming that the transmission rate is time dependent and specifically is determined by a simple sinusoidal function, this view has changed in recent years:

  • seasonally forced models of childhood infections now more often use a square wave

  • the transmission rate is assumed to be high during school terms and low at other times

\[\beta(t)=\beta_{0}(1+b_{1}Term(t))\]

where:

  • \(Term(t)\) is +1 during the school term and −1 at other times.

  • \(b_{1}\) to represent the amplitude of seasonality

Table 1: Timings of the major school holidays when Term = −1; during all other times Term = +1.
Holiday Model days Calender dates
Chirstmas 356-6 December 21 - January 6
Easter 100-115 April 10-25
Summer 200-251 July 19 - September 8
Autumn Half Term 300-307 October 27-November 3

We obtain 92 school holidays , 273 days of school => many more “+1” days than “−1” days => give rise to a mean transmission rate averaged over the year. To ensure R0 is constant irrespective of the precise forcing function used and the amplitude of seasonality. If there are \(D+\) days of school and \(D−\) holidays, then our forcing function would be:

\[\beta(t)=\frac{\beta_{0}}{\frac{1}{365}((1+b_{1})D_{+}+(1-b_{1})D_{-})}(1+b_{1}Term(t))\]

Code 
library(latex2exp)

beta_0 <- 1250
beta_1 <- 0.25
t <- seq(0,365,le = 1001)

## sinusoidal function
omega <- (2*pi)/365
beta_sin_func <- data.frame(x = t,
                            y = beta_0*(1 + beta_1*cos(omega*t)))

## term time tranmission function

schools_time <- c(0, 7, 100, 116, 200, 252, 300,308,356)
schools_open <- c(-1, 1, -1, 1,   -1,   1,   -1, 1,-1)

beta <- approx(
  schools_time,
  beta_0*(1 + beta_1*schools_open),
  xout = t,
  method = "constant",
  rule = 2)

## correct term time tranmission function
mean_beta <- 1/365*((1+beta_1)*273+(1-beta_1)*92)

beta_correct <- approx(
  schools_time,
  beta_0/mean_beta*(1 + beta_1*schools_open),
  xout = t,
  method = "constant",
  rule = 2)

SIR model with term-based forcing function 
sir_term <- odin2::odin({
  N <- S + I + R
  deriv(S) <- b*N - (beta_t*(I/N) + b)*S
  deriv(I) <-beta_t*(I/N)*S - (b + gamma)*I
  deriv(R) <- gamma*I - R*b

  # seasonality forcing
  beta_0 <- parameter(0.4)
  beta_1 <- parameter(0)
  
  school_open <- parameter()
  school_time <- parameter()
  school_data_dim <- parameter()
  dim(school_time) <- school_data_dim
  dim(school_open) <- school_data_dim
  
  term <- interpolate(school_time, school_open, "constant")
  beta_t <- beta_0*(1 + beta_1*term)
  
  # initialize starting population
  init_S <- parameter(9500)
  init_I <- parameter(500)
  gamma <- parameter(0.05)
  b <- parameter(0.05)
  
  initial(S) <- init_S
  initial(I) <- init_I
  initial(R) <- 0

})
helper function 
get_adjust_factor <- function(b1, formatted_schooldays){
  # compute number of school days for each term
  total_school_days <- sapply(formatted_schooldays$school_days, \(pair){
    pair$max - pair$min + 1
  })
  # compute number of holidays for each term
  total_holidays <- sapply(formatted_schooldays$holidays, \(pair){
    pair$max - pair$min + 1
  })
  
  adjust_factor <- (1/365)*((1+b1)*sum(total_school_days) + (1-b1)*sum(total_holidays))
  1/adjust_factor
}

format_schooldays <- function(school_time, school_open){
  prev <- 0
  school_days <- list()
  holidays <- list()
  sapply(1:length(school_time), \(i){
    if(school_open[i] == -1 && school_time[i]>0){
      school_days[[length(school_days) + 1]] <<- list(min = prev, max = school_time[i]-1)
      prev <<- school_time[i]
    }else if (school_open[i] == 1 && school_time[i]>0){
      holidays[[length(holidays) + 1]] <<- list(min = prev, max = school_time[i]-1)
      prev <<- school_time[i]
    }
  })
  
  # handle cases when the end of school time is not the end of the year
  if(prev < 365){
    if(school_open[length(school_open)] == -1){
      holidays[[length(holidays) + 1]] <- list(min = prev, max = 365)
    }else{
      school_days[[length(school_days) + 1]] <- list(min = prev, max = 365)
    }
  }
  
  list(
    school_days = school_days,
    holidays = holidays
  )
}

# get a list of item, each item is a list of 2, for min x and max x for the annotation
get_annotate_layers <- function(annotate_range, color = "red", alpha = 0.5){
  lapply(annotate_range, \(pair){
    annotate("rect", 
             xmin = pair$min, xmax = pair$max,
             ymin = -Inf, ymax = Inf,
             fill = color, alpha = alpha)
  })
}
unadjusted forcing function 
school_time <- c(0, 7, 100, 116, 200, 252, 300,308,356)
school_open <- c(-1, 1, -1, 1,   -1,   1,   -1, 1,-1)

term_pars <- list(
  school_time = school_time, 
  school_open = school_open,
  school_data_dim = length(school_time),
  beta_0 = 0.125,
  beta_1 = 0.29,
  gamma = 1/13,
  b = 0.02,
  init_S = 499900,
  init_I = 100
  )

# format school data for generating annotation layer
school_data <- format_schooldays(school_time, school_open)

# run model and plot output
run_mod(sir_term, term_pars, duration = 365, timestep = 1) %>% 
  ggplot() +
    geom_line(aes(x = t, y = I), color = "cornflowerblue") +
  get_annotate_layers(school_data$holiday, alpha = 0.3) +
  labs(
    title = "Number of infectives over a year",
    x = "Time (days)",
    y = "Number of infectives"
  )+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 18))

adjusted forcing function 
adjust_factor <- get_adjust_factor(term_pars$beta_1, school_data)

# use adjusted value for beta_0 instead
adjusted_pars <- term_pars
adjusted_pars$beta_0 <- adjust_factor*adjusted_pars$beta_0

# run model again and plot
run_mod(sir_term, adjusted_pars, duration = 365, timestep = 0.5) %>% 
  ggplot() +
    geom_line(aes(x = t, y = I), color = "cornflowerblue") +
  get_annotate_layers(school_data$holiday, alpha = 0.3) +
  labs(
    title = "Number of infectives over a year",
    x = "Time (days)",
    y = "Number of infectives"
  )+
  theme(axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        strip.text.x = element_text(size = 18),
        plot.title = element_blank(),
        plot.subtitle = element_text(size = 18),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 18))

Code materials : https://quingzz.github.io/blog/seasonality/seasonality.html

References

Bailey, Norman TJ. 1975. The Mathematical Theory of Infectious Diseases and Its Applications.
Becker, Alexander D., and Bryan T. Grenfell. 2017. “tsiR: An R Package for Time-Series Susceptible-Infected-Recovered Models of Epidemics.” PLOS ONE 12 (9): e0185528. https://doi.org/10.1371/journal.pone.0185528.
Silverman, Bernard W. 2018. Density Estimation for Statistics and Data Analysis. Routledge.